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A method is presented for the direct calculation of tunneling corrections for unsymmetrical Eckart type 
potential barriers. It is based on a modified 6-point Gaussian quadrature formula. Accuracy is better than 1 per- 
cent over a wide range of tunneling parameter values. 
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1. Introduction 

The Eckart potential function [l] 1 is often used to estimate 
quantum mechanical tunneling corrections to theoretically 
determined chemical rate constants. The correction factor 
r* is defined as the ratio of the quantum mechanical to the 
classical mechanical barrier crossing rate. It can be ex- 
pressed [2] as in integral over the energy E, 

T* = exp(VJkT) P Kexp(-E/kT)dE/kT (1) 

where V x is the height of the potential barrier, and K is the 
transmission probability for tunneling. K depends on E and 
three other parameters which are determined by the shape 
of the barrier and an effective mass for the system. Johnston 
and Heicklen [3] have evaluated this integral numerically 
for a number of parameter values. For certain applications 
their results are inconvenient to use because interpolation is 
required to get values not tabulated. In view of this, I have 
devised a simple method which can be used to calculate T* 
directly, for any set of parameter values within the ranges 
chosen by Johnston and Heicklen. The method is presented 
in the form of a small FORTRAN subroutine called TUNL. 
In the next section, the details of the method are discussed. 
Following this, the results of a series of comparisons with an 
accurate calculation are presented. Finally, the subroutine 
is listed in the Appendix. 



2. Derivation of the method 

Eckart's potential has the form 

V= -y[A-BI(l-yMl-y) 
y = — exp(27nr/Z,) 
A = V x -V 2 

b = (v»+ v$y 

l = 2*(-2iF'y*(rp+ v;*)- 1 

The potential has the limiting value of zero when x — — oo, 
goes through a single maximum of height V x as x increases, 
and has a limiting value of V x — V 2 as x — • + oo. F* is the 
second derivative of Fat its maximum. The lower bound E 
in the integral (1) is equal to zero when V x < V 2 , and to V x 
— V 2 when V x > V 2 . The three parameters used by 
Johnston and Heicklen are a lt a 2 , and u*. 



u* = hv'IkT 

a, = 2*VJhv*, i = 1,2 

v* = (U2icX-F'lmy* 



(2) 



where m is an effective mass for tunneling (see ref. 2, p. 53). 
The integral (1) can be written in a symmetrical form by in- 
troducing a new variable, e = (E — V x )lkT. It becomes 



T* = fKe^de 



(3) 
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where e„ = — v x = — VJkT when V x < V 2 , and e„ = — v 2 
= — VJkT when V x > V 2 . In terms of the parameters (2), 
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the transmission probability K, derived by Eckart, has the 
form 

K = [cosh27r(a! + a^ — cosh2ir(a l — aJ]X 
[cosh27r(ai + a 2 ) + DY 1 



2-ko, = 7r[(e + v,)lc]*, i = 1,2 
c = (l/SjTruV,-^ + a 2 "*) 2 
D = cosh2x</ if d is real 

= cos2x|</| if d is imaginary 
2zd = (4a,a 2 - tt 2 )* 



(4) 



The method used to evaluate (3) is a modified 6-point 
Gaussian quadrature formula based on Legendre polyno- 
mials [4]. This was used even though the nature of the in- 
tegral suggests using a formula based on Laguerre poly- 
nomials. When the number of evaluations of K is kept small, 
neither of these methods is satisfactory for the whole range 
of parameter values used by Johnston and Heicklen, so a 
variation of the first method was developed. 

When e gets large, K approaches unity. The method uses 
a Gaussian formula for that part of the integral where K < 
1. The remainder where K ~ 1 is evaluated analytically. 
Thus, if K{e) = 1 for e > e 6 , then 

\~ K(t)e-*de « (°Vde = e-\ 

To evaluate e b , examine (4) as e — oo. One gets a { — 
V2{elcy* and 

K - 1 - (1 + Z?)(i/ 2 exp(27r(e 6 /c)*) + #)" = K b . 

Setting K b to some value close to unity and solving this 
equation for e b gives 



— firMf^])- 



(5) 



It happens that this value is not entirely satisfactory, and 
subtracting from it the average value of v x and v 2 gives bet- 
ter results. Also, in some cases, e b calculated in this way is 
very large. There is no point in using this value for e* as the 
upper bound of the Gaussian formula if the integrand at 
this point is negligible because of the exponential factor. 
Thus e b was kept below a certain fixed value e m „. There 
resulted two parameters, K b and e m „, which were adjusted 
to minimize the sum of the squares of the differences be- 
tween the results of this method and the corresponding 
tabulated values of Johnston and Heicklen. 



3. Test of the method 

Extensive testing of the accuracy of the method was per- 
formed by comparing it with an accurate 40-point Gaussian 
formula having the cutoff fixed at e b corresponding to K b = 
0.999 or at e b < 8. In the ranges 0.5 < c^ < a 2 < 20, and 
2 < u* < 16, a group of 10,910 comparisons was made. For 
this set there were the additional restrictions that when a! 
> 8 then u* < 12, or when a x > 16 then u* < 10. Note 
that r*(a n aj = r*(a 2 , a x ). A second set of 4,920 compar- 
isons was made in the ranges 0.5 < c^ < a 2 < 20, and 0.05 
< u* < 1.5. The results of these tests are given in table 1 
in the form of histograms. These show the number of values 
which differ from the accurate values by a given percentage 
range. It can be seen from these results that very few values 
are in error by as much as 5 percent. Such accuracy should 
be quite adequate for most rate constant calculations. 



TABLE I. Tests of the accuracy of TUNL 



Variation 


from 


Number of differences in 


accurate i 


values 


percentage ranges 


Percent difference 


Set I - * 


Set II - c 


-5.5%, -4.5% 




3 


1 


-4.5, -3.5 




4 


5 


-3.5, -2.5 




6 


19 


-2.5, -1.5 




30 


26 


-1.5, -0.5 




3246 


140 


-0.5, 0.5 




6811 


4475 


0.5, 1.5 




343 


227 


1.5, 2.5 




217 


27 


2.5, 3.5 




170 




3.5, 4.5 




78 




4.5, 5.5 




2 




Standard deviations 


0.77% 


0.42% 



* For both Sets I and II, a 2 > or, . The values used were 0.5, 1.0, 1.5, 2.0, 
. . . 40.0. 

* For Set I, u* = 2, 3,4, . . . 16. Also if ct x >8thenu*< 12 and if a, > 
16 then u* < 10. 

' For Set II, u* = 0.05, 0.1, 0.2, 0.5, 1.0, 1.5. 



4. Appendix. Listing of TUNL 

The parameters for this program are ALPHl = a x , 
ALPH2 = a 2 , U = u* t and G = r*. It will calculate T* ac- 
curately in the parameter ranges 0.5 < o^ , ct 2 < 20, and 
< u* < 16 with the additional restrictions that when c^ 
and a 2 > 8 then u* < 12, and when a x and cc 2 > 16 then 
u* < 10. 
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SDBROUTINE TITNL ( ALPHI . ALPH2 , V. G) 

DIMENSION X(6) ,W(6) 

DATA X/-. 9324695. -.6612094. -.2386192.. 2386192, 

DATA W/. 1713245.. 3607616,. 4679139.. 4679139,. 3607* 

DATA PI, PI 2, PI SQ/3. 1415926, 6. 2831853, 9. 8696044/ 

C0SH(Z)=.5»(EXP(Z)+EXP(-Z) ) 

DPI2=U/PI2 

C=. 1 2 5 •PI«U«(1./SQRT( ALPHI ) +1 . / SQRT ( ALPH2 ) ) ««2 

V1-=UPI2«ALPH1 

V2-UPI2»ALPH2 

D=4 . «ALPH1 «ALPH2-PISQ 

IF(D.LT.O) GOTO 10 

DF=COSH(SQRT(D)> 

GOTO 11 

10 DF=COS(SQRT(-D)) 

11 IF(V2.GE.V1) EZ — VI 
IF(V1.GT.V2) EZ=-V2 

EB-AMINl(C*(ALOG(2.*(l.+DF)/.014)/PI2)**2-.5*(Vl<i 
EM=.5«(EB-EZ) 

EP= .5MEB+EZ) 
G=0 

DO 20 N=l ,6 
E-EM«X(N)+EP 
A1 = PI«SQRT( (E+VD/C) 
A2=PI«SQRT( (E+V2)/C) 
FP=C0SH(A1+A2) 
FM=C0SH(A1-A2) 
20 G=G+W(N)»EXP(-E)«(FP-FM)/(FP+DF) 
G=EM«G+EXP(-EB) 
RETURN 
END 
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